/*
 * Qmc.cpp
 *
 *  Created on: Jan 24, 2010
 *      Author: snirgaz
 */

#include "Qmc.h"

void Qmc::init(QmcParameters &qmcParams) {
	qmcParams_ = qmcParams;
#ifdef PT
	String OutFile = qmcParams.outPath;
	OutFile.append("/Conf").append(
			boost::lexical_cast<string>(qmcParams.confNum)).append(".out");
#else
	string OutFile = qmcParams.outPath.getVal();
	OutFile.append("/Conf").append(".out");
#endif
	outPutFile_.open(OutFile.c_str(), ios::out);
	outPutFile_ << "Num Of Cores :" << qmcParams.numOfCores.getVal() << endl;
	outPutFile_ << "Seed Number :" << qmcParams.seed.getVal() << endl;
#ifdef FORCE
	qmcParams_.greenFactorV = (2. / 3.) * qmcParams.epsilon;
	qmcParams_.greenFactorF = (1. / 9.) * pow(qmcParams.epsilon, 3)
	* (2. * qmcParams.lambda);
#else
	qmcParams_.greenFactorV.setVal(qmcParams.epsilon.getVal());
	qmcParams_.greenFactorF.setVal(0.);
#endif
//	qmcParams.sigmaJumpAccept = sqrt(
//			qmcParams.jumpAcceptFactor * qmcParams.mbar) * qmcParams.sigma;
#ifdef NOSWAP
	qmcParams_.swapProb = 0;
#endif
	qmcParams_.mbarOpen = qmcParams_.mbarSwap = qmcParams_.mbarWiggle =
			qmcParams_.mbar;
	// Simluation params

}
void Qmc::initConf() {
	if (qmcParams_.init.getVal())
		this->initializeFromFile();
	else {
#ifdef POT_AZIZ
		string init_type = qmcParams_.initialType.getVal();
//		if (init_type == "VACANCIES")
//			wl_.initializeSolidBCCVacancies(qmcParams.numOfVacancies,
//					qmcParams.outPath, static_cast<QmcRnd&>(*this));
//		if (init_type == "INTERSTITALS")
//			wl_.initializeSolidBCCInterstitials(
//					qmcParams.numOfInterstitals, qmcParams.outPath,
//					static_cast<QmcRnd&>(*this));
//		if (init_type == "HCP")
//			wl_.initializeSolidHCP(qmcParams.outPath,
//					qmcParams.latFileInit);
//		if (init_type == "GRAIN_B") {
//			wl_.initializeSolidBCCFixedGB(qmcParams.outPath,
//					qmcParams.latFileInit, qmcParams.NumPFixed);
//			qmcParams_.NumPMoving.setVal(qmcParams.NumP - qmcParams.NumPFixed);
//			QmcRnd::particleNum_.updateParams(
//					UniIntParams(0, qmcParams.NumPMoving - 1));
//			String fileName = qmcParams.outPath + "/numOfFixed.xml";
//			pugi::xml_document doc;
//			pugi::xml_node root_node = doc.append_child("OUTPUT_PARAMETERS");
//			addData(root_node, "NUM_FIXED", qmcParams.NumPFixed);
//			doc.save_file(fileName.c_str());
//		}
//		if (init_type == "REGULAR_BCC")
//			WorldLines::initializeSolidBCC(qmcParams.outPath);
//		if (init_type == "REGULAR")
//			WorldLines::initializeSolid(qmcParams.outPath);
//#endif
//#ifdef POT_LOG
//		WorldLines::initializeSolid(qmcParams.outPath);
//#endif

		energyCalc_.calcConfigurationEnergy();
//		this->printEnergyDiff();
	}
}

void Qmc::initializeWLDataFromFile(String fileName) {
	wl_.readFromFile(fileName);
}
void Qmc::initializeFromFile() {
	string dir(qmcParams_.initPath.getVal()), fileName;
	outPutFile_ << "Initialize From File";
	// Read  WL data
	Qmc::initializeWLDataFromFile(dir + "/WL.h5");
	// Read previous sim xml File
	fileName = dir + "/simParams.xml";
	pugi::xml_document xmlDoc;
	pugi::xml_node rootNode, simParam;
	try {
		pugi::xml_parse_result result = xmlDoc.load_file(fileName.c_str());
		if (!result)
			throw result;
	} catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	rootNode = xmlDoc.child("SIM_PARAMETERS");
	qmcParams.gamma = readXmlData<double>(rootNode, "GAMMA");
	qmcParams.rate = readXmlData<int>(rootNode, "SWEEP_RATE");
	qmcParams.mbarSwap = readXmlData<int>(rootNode, "MBAR_SWAP");
	qmcParams.mbarOpen = readXmlData<int>(rootNode, "MBAR_OPEN");
	qmcParams.mbarWiggle = readXmlData<int>(rootNode, "MBAR_WIGGLE");
	// Load stream Data
	fileName = dir + "/randomStream";
	QmcRnd::loadStream(fileName);
	// Print initial configuration
	String OutFile = qmcParams.outPath;
	OutFile.append("/WLStart");
	WorldLines::dumpToFile(OutFile);
	Measure::readMeasurements(dir + "/Measure.h5", outPutFile_);
}
void Qmc::thermalize() {
	outPutFile_ << "Configuration Number - " << qmcParams.confNum << endl;
	timeval startTime, currentTime, result;
// Start Time
	gettimeofday(&startTime, NULL);
// Thermalize
	outPutFile_ << "Thermalization of " << qmcParams.numThermal << " :  \n"
			<< std::endl;
// Adjust gamma
	this->adjustGamma();
// Check energy difference
	this->printEnergyDiff();
	gettimeofday(&currentTime, NULL);
	timersub(&currentTime, &startTime, &result);
	outPutFile_ << "\n Thermalization time " << (double(result.tv_sec)) << "\n"
			<< std::endl;
}
void Qmc::runN(int n) {
	ippSetNumThreads(1);
	for (qmcParams.step = 1; (qmcParams.step < n); qmcParams.step++) {
		isAccepted_ = this->nextMove();
	}
}
bool Qmc::run() {
	TimePoint startTime = boost::chrono::system_clock::now();
	Hour dumpTime(qmcParams.dumpTime);
	Hour finishTime(qmcParams.finishTime);
	bool done = true;
	outPutFile_ << "Dump Every " << dumpTime.count() << " Hour" << std::endl;
	outPutFile_ << "Finish After " << finishTime.count() << " Hour"
			<< std::endl;
	this->dump_params();
	this->initialize();
	if (!(qmcParams.init))
		this->thermalize();
	outPutFile_ << "Sweep Rate" << qmcParams.rate << std::endl;
// Initialize WL
	this->initializeCountOutput();
// Start Time
	for (qmcParams.step = 1; (!Measure::stopMeasure()); qmcParams.step++) {
		isAccepted_ = this->nextMove();
		//	if (is_accepted_) this->dump(qmcParams.dirnum);
		bool flagSweep = !(qmcParams.step % qmcParams.rate);
		if (flagSweep)
			Measure::performMeasure();
		// Dump after 1 Hour
		Hour timeDiff = boost::chrono::system_clock::now() - startTime;
		if (timeDiff.count() > dumpTime.count()) {
			this->dump();
			dumpTime += Hour(qmcParams.dumpTime);
		}
		// Stop Simulation
		if (timeDiff.count() > finishTime.count()) {
			done = false;
			break;
		}
	}
	Hour timeDiff = boost::chrono::system_clock::now() - startTime;
	outPutFile_ << "\n Measurements time " << timeDiff.count() << "\n"
			<< std::endl;
// Check energy difference
	this->printEnergyDiff();
	return done;
}
bool Qmc::nextMove() {
//	assert(this->loop_check(qmcParams.step));
//	assert(WL_->isValidConf(qmcParams.outpath, B_->get_walk(), SP_));
//	static int i = 0;
//	if (qmcParams.flag) {
//		String file_name = qmcParams.outPath + "/WL"
//				+ boost::lexical_cast<String>(i) + ".h5";
//		WorldLines::dumpToFile(file_name);
//		file_name = qmcParams.outPath + "/RW" + boost::lexical_cast<String>(i)
//				+ ".h5";
//		MoveParams::dumpBeads(file_name);
//		i++;
//	}
	assert(this->EKcheck(qmcParameters_.step, moveParams_));
	assert(this->energy_diff(qmcParameters_.step));
	int next;
	if (wl_.wlData_.getOpenClose() == CLOSED)
		next = 0;
	else {
		// draw next move
		next = nextMove_.draw();
	}
	moveParams_->moveType = MoveType(next);
	stats_.moves[next];
// Prepare next move
	switch (moveParams_->nextMove_) {
	case OPEN_MOVE: //WIGGLE
		return openMove_.makeMove();
	case CLOSE_MOVE: //CLOSE
		return closeMove_.makeMove();
	case WIGGLE_MOVE: //SWAP
		return wiggleMove_.makeMove();
	case SWAP_MOVE: //OPEN
		return swapMove_.makeMove();
	default:
		return false;
	}

}

double Qmc::calcRhoOpenClose() {
	double sigmaMove = qmcParams.sigma * sqrt(double(MoveParams::mbar - 1)), C =
			-0.5 * pow(sigmaMove, -2);
	BeadPos start, end, delta;
	start = WL::getDataPos(MoveParams::wlStart, MoveParams::mStart);
	end = WL::getDataPos(MoveParams::wlEnd, MoveParams::mEnd);
	end.subBeadPos(delta, start);
	double drSqr = delta.sumOfSqures();
	return exp(C * drSqr);
}

void Qmc::dump() {
	string dir(qmcParams.outPath), fileName, param;
// Dump World Lines
	fileName = dir;
	fileName.append("/WL.h5");
	WorldLines::dumpToFile(fileName);
	fileName = dir;
// Dump Beads
	fileName.append("/RW.h5");
	MoveParams::dumpBeads(fileName);
// Dump Stats
	fileName = dir;
	fileName.append("/Stats.xml");
	stats_.dumpToFile(fileName);
// Save stream position
	fileName = dir + "/randomStream";
	QmcRnd::saveStream(fileName);
// Dump Measurements.
	fileName = dir;
	fileName.append("/Measure.h5");
	Measure::dumpMeasurements(fileName, outPutFile_);
// Write Parameters to XML file
	fileName = dir;
	fileName.append("/Output_Params.xml");
	pugi::xml_document doc;
	pugi::xml_node root_node = doc.append_child("OUTPUT_PARAMETERS");
	addData(root_node, "NEXT_MOVE", MoveParams::next);
	addData(root_node, "WORM_WL_START", MoveParams::wlStart);
	addData(root_node, "WORM_WL_END", MoveParams::wlEnd);
	addData(root_node, "WORM_WL_SWAP_START", MoveParams::wlSwapedStart);
	addData(root_node, "WORM_WL_SWAP_END", MoveParams::wlSwapedEnd);
	addData(root_node, "M_BAR", MoveParams::mbar);
	addData(root_node, "MSTART", MoveParams::mStart);
	addData(root_node, "MEND", MoveParams::mEnd);
	addData(root_node, "TAU", MoveParams::tau);
	addData(root_node, "IS_CLOSED", WL::getOpenClose() == CLOSED);
	doc.save_file(fileName.c_str());
// Parameters for next simulation
	fileName = dir;
	fileName.append("/simParams.xml");
	pugi::xml_document doc2;
	root_node = doc2.append_child("SIM_PARAMETERS");
	addData(root_node, "NUM_CLOSED_MEASURE", Measure::getNumClosedMeasure());
	addData(root_node, "NUM_OPEN_MEASURE", Measure::getNumOpenMeasure());
	addData(root_node, "GAMMA", qmcParams.gamma);
	addData(root_node, "SWEEP_RATE", qmcParams.rate);
	addData(root_node, "MBAR_SWAP", qmcParams.mbarSwap);
	addData(root_node, "MBAR_OPEN", qmcParams.mbarOpen);
	addData(root_node, "MBAR_WIGGLE", qmcParams.mbarWiggle);
	addData(root_node, "DONE", Measure::stopMeasure() ? 1 : 0);
	doc2.save_file(fileName.c_str());
}
void Qmc::dump_params() {
	string fileName(qmcParams.outPath), param;
	fileName.append("/SP_params.xml");
	pugi::xml_document doc;
	pugi::xml_node root_node = doc.append_child("SIM_PARAMETERS");
	addData(root_node, "EPSILON", qmcParams.epsilon);
	addData(root_node, "SIGMA", qmcParams.sigma);
	doc.save_file(fileName.c_str());
}
void Qmc::updatePath(String &path) {
	qmcParams.outPath = path;
}
void Qmc::updateRS(double rs) {
#ifdef ROT
	qmcParams.lambda = 1
	/ (PI * double(qmcParams.NumP) * pow(rs, 2) / (sqrt(3) / 2.));
#else
	qmcParams.lambda = 1 / (PI * double(qmcParams.NumP) * pow(rs, 2));
#endif
	this->updateLambda(qmcParams.lambda);
	WorldLines::EnergyData::calcEnergy(qmcParams.greenFactorV,
			qmcParams.greenFactorF);
}
void Qmc::updateLambda(double lambda) {
	qmcParams.lambda = lambda;
	qmcParams.sigma = sqrt(2 * qmcParams.lambda * qmcParams.epsilon);
#ifdef FORCE
	qmcParams.greenFactorV = (2. / 3.) * qmcParams.epsilon;
	qmcParams.greenFactorF = (1. / 9.) * pow(qmcParams.epsilon, 3)
	* (2. * qmcParams.lambda);
#else
	qmcParams.greenFactorV = qmcParams.epsilon;
#endif
	QmcRnd::normal_.updateParams(NormalVariableParams(0, qmcParams.sigma));
}
void Qmc::adjustGamma() {
	String fileName;
	outPutFile_ << "Adjust Gamma: \n" << std::endl;
	qmcParams.mbarSwap = qmcParams.mbarOpen = qmcParams.mbarWiggle =
			qmcParams.mbar;
	int max_mbar = qmcParams.mbar + 2, min_mbar = 4; // max(4, qmcParams.mbar / 2);
	DoubleVec swaps, opens, wiggles, closes, gammas(5), open_prob(5);
	boost::scoped_ptr<boost::progress_display> gamma_progress(
			new boost::progress_display(qmcParams.numThermal, outPutFile_));
	int moves_per_sweep = qmcParams.numThermal / 5;
	qmcParams.mbarSwap = qmcParams.mbarOpen = qmcParams.mbarWiggle = min_mbar;
	QmcRnd::posWiggle_.updateParams(
			UniSmallIntParams(1, qmcParams.mbarWiggle - 2));
	QmcRnd::posOpen_.updateParams(UniSmallIntParams(1, qmcParams.mbarOpen - 2));
	for (int j = 0; j < 5; j++) {
		for (qmcParams.step = 1; qmcParams.step < moves_per_sweep;
				qmcParams.step++, ++(*gamma_progress)) {
			if ((qmcParams.step) >= 5)
				qmcParams.flag = true;
			else
				qmcParams.flag = false;
			qmcParams.flag = false;
			isAccepted_ = this->nextMove();
		}
		int total_prob = 0;
		for (int i = 0; i < 4; i++)
			total_prob += stats_.moves[i];
		open_prob[j] = double(stats_.moves[3]) / double(total_prob);
		gammas[j] = qmcParams.gamma;
		qmcParams.gamma /= open_prob[j] / qmcParams.openProb;
		stats_.initalize();
	}
	for (int j = 0; j < 5; j++)
		outPutFile_ << "\n Sweep :" << j << " Gamma : " << gammas[j]
				<< " Open_prob : " << open_prob[j] << std::endl;
	gamma_progress.reset(
			new boost::progress_display(qmcParams.numThermal, outPutFile_));
	for (int j = 0; j < 5; j++) {
		for (qmcParams.step = 1; qmcParams.step < moves_per_sweep;
				qmcParams.step++, ++(*(gamma_progress.get()))) {
			isAccepted_ = this->nextMove();
		}
		int total_prob = 0;
		for (int i = 0; i < 4; i++)
			total_prob += stats_.moves[i];
		open_prob[j] = double(stats_.moves[3]) / double(total_prob);
		gammas[j] = qmcParams.gamma;
		qmcParams.gamma /= open_prob[j] / qmcParams.openProb;
		stats_.initalize();
	}
	fileName = qmcParams.outPath + "/WLAdjust.h5";
	WorldLines::dumpToFile(fileName);
	for (int j = 0; j < 5; j++)
		outPutFile_ << "\n Sweep :" << j << " Gamma : " << gammas[j]
				<< " Open_prob : " << open_prob[j] << std::endl;
	int mbar = min_mbar;
	while (mbar <= max_mbar) {
		qmcParams.mbarSwap = qmcParams.mbarOpen = qmcParams.mbarWiggle = mbar;
		QmcRnd::posWiggle_.updateParams(
				UniSmallIntParams(1, qmcParams.mbarWiggle - 2));
		QmcRnd::posOpen_.updateParams(
				UniSmallIntParams(1, qmcParams.mbarOpen - 2));
		for (qmcParams.step = 1; qmcParams.step < moves_per_sweep;
				qmcParams.step++)
			isAccepted_ = this->nextMove();
		// Adjust gamma
		int total_prob = 0;
		for (int i = 0; i < 4; i++)
			total_prob += stats_.moves[i];
		open_prob[0] = double(stats_.moves[3]) / double(total_prob);
		qmcParams.gamma /= open_prob[0] / qmcParams.openProb;
		stats_.initalize();
		for (qmcParams.step = 1; qmcParams.step < moves_per_sweep;
				qmcParams.step++)
			isAccepted_ = this->nextMove();
		swaps.push_back(
				fabs(
						0.5
								- double(stats_.acceptedBoltzman[2])
										/ double(stats_.moves[2])));
		opens.push_back(
				fabs(
						0.5
								- double(stats_.acceptedBoltzman[3])
										/ double(stats_.moves[3])));
		wiggles.push_back(
				fabs(
						0.5
								- double(stats_.acceptedBoltzman[0])
										/ double(stats_.moves[0])));
		closes.push_back(
				fabs(
						0.5
								- double(stats_.acceptedBoltzman[1])
										/ double(stats_.moves[1])));
		outPutFile_ << "mbar = " << mbar << "swaps = " << swaps.back()
				<< " opens = " << opens.back() << " wiggle = " << wiggles.back()
				<< " closed = " << closes.back() << std::endl;
		mbar++;
	}
	int delta_m;
	delta_m = min_element(swaps.begin(), swaps.end()) - swaps.begin();
	qmcParams.mbarSwap = min_mbar + delta_m;
	delta_m = min_element(opens.begin(), opens.end()) - opens.begin();
	qmcParams.mbarOpen = min_mbar + delta_m;
	delta_m = min_element(wiggles.begin(), wiggles.end()) - wiggles.begin();
	qmcParams.mbarWiggle = min_mbar + delta_m;
	QmcRnd::posWiggle_.updateParams(
			UniSmallIntParams(1, qmcParams.mbarWiggle - 2));
	QmcRnd::posOpen_.updateParams(UniSmallIntParams(1, qmcParams.mbarOpen - 2));
	outPutFile_ << "Final mbar Swap : " << qmcParams.mbarSwap << std::endl;
	outPutFile_ << "Final mbar Open : " << qmcParams.mbarOpen << std::endl;
	outPutFile_ << "Final mbar Wiggle : " << qmcParams.mbarWiggle << std::endl;
	for (int j = 0; j < 5; j++) {
		for (qmcParams.step = 1; qmcParams.step < moves_per_sweep;
				qmcParams.step++)
			isAccepted_ = this->nextMove();
		// Adjust gamma
		int total_prob = 0;
		for (int i = 0; i < 4; i++)
			total_prob += stats_.moves[i];
		open_prob[j] = double(stats_.moves[3]) / double(total_prob);
		gammas[j] = qmcParams.gamma;
		qmcParams.gamma /= open_prob[j] / qmcParams.openProb;
		stats_.initalize();
	}
	for (int j = 0; j < 5; j++)
		outPutFile_ << "\n Sweep :" << j << " Gamma : " << gammas[j]
				<< " Open_prob : " << open_prob[j] << std::endl;
}
void Qmc::reCalcEnergy() {
	WorldLines::calcEnergy(qmcParams.greenFactorV, qmcParams.greenFactorF);
}

void Qmc::updateRs(double rs) {
#ifdef ROT
	double lambda = 1
	/ (PI * double(qmcParams.NumP) * pow(rs, 2) / (sqrt(3) / 2.));
#else
	double lambda = 1 / (PI * double(qmcParams.NumP) * pow(rs, 2));
#endif
	this->updateLambda(lambda);
}
bool Qmc::energy_diff(int step) {
	EnergyData ED1, ED2;
	ED1 = static_cast<EnergyData&>(WorldLines::getEnergyData());
	ED2 = EnergyCalc::calcConfigurationEnergy();
	double diff = fabs(ED1.getEnergy() - ED2.getEnergy());
//ED2 = WorldLines::getEnergyData();
	if (diff > 0.001) {
		cout << step;
		return false;
	}
	return true;
}
void Qmc::printEnergyDiff() {
	EnergyData ED1, ED2;
	ED1 = static_cast<EnergyData&>(WorldLines::getEnergyData());
	ED2 = EnergyCalc::calcConfigurationEnergy(false);
	outPutFile_ << "\n Energy Difference " << ED1.getEnergy() - ED2.getEnergy()
			<< "\n" << std::endl;
	outPutFile_ << "\n Energy " << ED1.getEnergy() << "   " << ED2.getEnergy()
			<< "\n" << std::endl;
	outPutFile_ << "\n Energy V Difference "
			<< ED1.getEnergyV() - ED2.getEnergyV() << "\n" << std::endl;
	outPutFile_ << "\n Energy V " << ED1.getEnergyV() << "   "
			<< ED2.getEnergyV() << "\n" << std::endl;
	outPutFile_ << "\n Energy F Difference "
			<< ED1.getEnergyF() - ED2.getEnergyF() << "\n" << std::endl;
	outPutFile_ << "\n Energy F " << ED1.getEnergyF() << "   "
			<< ED2.getEnergyF() << "\n" << std::endl;
	double diff = ED1.getEnergy() - ED2.getEnergy();
}

bool Qmc::EKcheck(int step, int move) {
	double ek1, ek2;
	ek1 = WorldLines::getEnergyK();
	ek2 = WL::calcEK();
	if (fabs(ek1 - ek2) > 0.00001) {
		outPutFile_ << "\n";
		outPutFile_ << "From Data : " << ek1 << "\n";
		outPutFile_ << "From Calc : " << ek2 << "\n";
		outPutFile_ << "Step : " << step << "\n";
		outPutFile_ << "Move : " << move << "\n";
		outPutFile_ << std::endl;
		return false;
	} else
		return true;
}
void Qmc::swapWL(Qmc &Qmc2) {
	WorldLines &WL1 = static_cast<WorldLines&>(*this);
	WorldLines &WL2 = static_cast<WorldLines&>(Qmc2);
	WL1.swapData(WL2);
	this->reCalcEnergy();
	Qmc2.reCalcEnergy();
}
void Qmc::getEnergyCalcEnergy(WLEnergyData &E1, WLEnergyData &E2) {
// Current
	E1 = WorldLines::getEnergyData();
	EnergyData EP = EnergyCalc::calcConfigurationEnergy(false);
// Calc
// Kinetic Energy
	E2.setEnergyK(WL::calcEK());
	E2.setEnergyP(EP);
}
bool Qmc::measureDone() {
	return Measure::stopMeasure();
}
void Qmc::setNumOfThermal(int numThermal) {
	qmcParams.numThermal = numThermal;

}
bool Qmc::isClosed() {
	return WorldLines::isOpenClosed_ == CLOSED;
}
int Qmc::getSweepRate() {
	return qmcParams.rate;
}
double Qmc::getEK() {
	return WorldLines::getEnergyK();
}
WLEnergyData Qmc::getED() {
	return WorldLines::getEnergyData();
}
void Qmc::performMeasure() {
	Measure::performMeasure();
}
void Qmc::initializeCountOutput() {
	outPutFile_ << "Measurements :" << endl;
	Measure::initializeCountOutput(outPutFile_);
}
